####################################################################
### Replication Data  for:                                         #
### Hemesath S., Tepe M. (2024) Multidimensional preferences for   #
### technology risk regulation. The role of political beliefs,     #
### technology attitudes, and national innovation cultures".       #
### Regulation & Governance. DOI: xxxx                             #
####################################################################

### Required Packages

library(tidyverse)
library(Hmisc)
library(qualtRics)
library(reshape2)
library(foreign)
library(ggmap)
library(ggeffects)
library(cjoint)
library(cregg)
library(psych)
library(corrplot)
library(PerformanceAnalytics)
library(ggrepel)
library(ggridges)
library(broom)
library(ggh4x)
library(cowplot)
library(survey)

rm(list=ls())

########################
### Load Conjoint Data##
########################

CJ <- readRDS("Hemesath_Tepe_RegGov_2024_Replication_AVConjoint.RDS")

##################################
## Build and test index measures##
##################################

#############################
## I. Technology Attitudes ##
#############################

cols <- c("tech1","tech2","tech3","tech4","tech5",
          "tech6","tech7","tech8","tech9","ADMmakeeasy","ADMpreferhuman","ADMuneasy","ADMlessdiscr") # select items
CJ[cols] <- sapply(CJ[cols],as.numeric) # transform to numeric

#Recode reverse coded items
CJ$tech2 <- 6-CJ$tech2
CJ$tech3 <- 6-CJ$tech3
CJ$tech5 <- 6-CJ$tech5
CJ$ADMpreferhuman <- 6- CJ$ADMpreferhuman
CJ$ADMuneasy <- 6- CJ$ADMuneasy

#subset tech items per respondent
Tech <- select(CJ,all_of(cols))

names(Tech) <- c("tech1","tech2","tech3","tech4","tech5",
                 "tech6","tech7","tech8","tech9","ADM1","ADM2","ADM3","ADM4")

# test for factorability
KMO(Tech) # Factorability
det(cor(Tech))
cor.plot(Tech) # tech2 and tech3 seem to not fit well
fa.parallel(Tech,fa="fa") # scree plot suggests four to six factors
nfactors(Tech) # model with six factors leads to rejection of Chisq test

#run FA with 6 factors
fatech <- fa(Tech,nfactors=6,n.obs=nrow(Tech),rotate="promax",max.iter=5000,SMC=T,fm="ml")

#plot results of factor analysis: Appendix Figure 2
pdf("Appendix_Figure2.PDF",width=10,height=6)
fa.diagram(fatech) # results are theoretically plausible
dev.off()

fa.diagram(fatech)

# determine consistency of additive scale based on items loading on correlating factors
# (identified as "Technology Optimism")

alpha(Tech[,-c(2:3,5,11,12,13),]) # acceptable internal consistency

# build index
#Build Index from 6 Items
CJ$techatt <- as.numeric((rowSums(CJ[,c("ADMmakeeasy","tech1","tech6","tech4","tech7","tech8","tech9"),])/7))

###############################
## Recode Covariate meausures##
###############################

# Bin continuous measures into three terciles

CJ$ecoLR <- as.numeric(CJ$ecoLR)
CJ$age <- as.numeric(CJ$age)
CJ$climateeco <- as.numeric(CJ$climateeco)
CJ$culLR <- as.numeric(CJ$culLR)
CJ$ADMuneasy <- as.numeric(CJ$ADMuneasy)


CJ <- CJ %>% mutate(
  techatt_bin = factor(cut_number(techatt,3),labels = c("Low","Medium","High")),
  climateeco_bin = factor(cut_number(climateeco,3),labels = c("Low","Medium","High")),
  ecoLR_bin = factor(cut_number(ecoLR,3),labels = c("Low","Medium","High")),
  culLR_bin = factor(cut_number(culLR,3),labels = c("Low","Medium","High")),
  age_bin = factor(cut_number(age,3),labels = c("Young","Medium","Old"))
         )

###################################
######## Empirical Analysis #######
###################################

# store CJ attributes in vector for convenience
att <- c("Safety.Standards","Safety.Approval","Moral","Legal","Data","Supervision")

# set Reference Levels for estimation of AMCE

CJ <- CJ %>% mutate(
  Safety.Approval = fct_relevel(Safety.Approval,"Self",after=0),
  Safety.Standards = fct_relevel(Safety.Standards,"identical",after=0),
  Moral = fct_relevel(Moral,"occupant",after=0),
  Legal = fct_relevel(Legal,"Owner",after=0),
  Data = fct_relevel(Data,"Driver",after=0),
  Supervision = fct_relevel(Supervision,"none",after=0)
)

# Create vector for renamed attributes and levels

ats <- c("Safety Standards\nvs. current",
         "Safety Approval\nby",
         "Ethical\nPrioritization",
         "Liability for\nAccidents",
         "Access to\nTelemtry Data",
         "Human Supervision\nof Autopilot")

lvl <- c("Identical","Much more strict","Stricter",
         "Self-Regulation","Public Authority","Independent Body",
         "Protect Occupants","Choice of Driver","Protect Others",
         "Owner","Manufacturer",
         "Only Driver","Public Authorities","Companies",
         "Not required","Constant","After Warning")

lvls <- list("Safety.Approval" = levels(CJ$Safety.Approval),
             "Safety.Standards" = levels(CJ$Safety.Standards),
             "Moral" = levels(CJ$Moral),
             "Legal" = levels(CJ$Legal),
             "Data" = levels(CJ$Data),
             "Supervision" = levels(CJ$Supervision)
             )

##########################################
## Figure 1: Estimate Baseline Results ###
##########################################

## Pooled Regressions
m1 <- setNames(lapply(att,function(x){
  svyglm(as.formula(paste("selected~",x)),design=svydesign(data=CJ,id=~ResponseId))
}),att)

## Pooled Marginal Means
mm1 <- bind_rows(setNames(lapply(att,function(x){
  ggeffect(m1[[x]],paste(x))
}),att),.id="feature")

## Pooled AMCEs
amce1 <- bind_rows(setNames(lapply(att,function(x){
  summary(margins::margins(m1[[x]]))
}),att),.id="feature")

##Create Dummy values for reference categories
amce1$factor <- str_replace(amce1$factor,amce1$feature,"")
b <- bind_rows(lapply(att,function(x){
  expand.grid(colnames(CJ[paste(x)]), levels(CJ[[paste(x)]]))
}))
names(b) <- c("feature","factor")
amce1 <- select(merge(b, amce1, by="factor", all = T),-c("feature.y"))
amce1[is.na(amce1)] <- 0

#prepare for plotting

dfa <- select(as.data.frame(mm1),c(1:3,5:6))
dfb <- select(amce1,c(1:3,7:8))
names(dfb) <- c("level","feature","predicted","lower","upper")
names(dfa) <- c("feature","level","predicted","lower","upper")
df <- bind_rows(dfa,dfb,.id="Analysis") %>%
  mutate(.,feature=factor(feature,levels=att,labels=ats),
         level=factor(level,levels=levels(level),labels=lvl),
         Analysis = factor(Analysis,levels=c(1,2),labels=c("Marginal Means",
                                                           "Average Marginal Component Effect")))
df$Analysis <- fct_relevel(df$Analysis,"Marginal Means",after=0)
df$level <- fct_relevel(df$level,"Much more strict",after=20)
levels(df$level)
##Plot Pooled Baseline Results

library(showtext)
font_families_google()
font_add_google("Roboto Condensed", family = "Roboto Condensed")
showtext_auto()

fig1 <- ggplot(df, aes(x = level, y = predicted))+
  geom_hline(aes(yintercept=ifelse(Analysis=="Marginal Means",0.5,0)),size=1,color="grey")+
  geom_pointrange(aes(ymin=lower,
                      ymax=upper),
                  position=position_dodge2(width=.75),
                  fatten=2.5,
                  size=1,
                  color="navy")+
  facet_grid(feature~Analysis,
             scales="free",
             switch="y",
             space = "free")+
  scale_x_discrete(position = "top",
                   expand=c(-.65,0))+
  coord_flip()+
  theme_classic(base_family = "sans")+
  theme(text=element_text(size=16,family="Roboto Condensed"),
        strip.text.y.left=element_text(angle=0,hjust=0,
                                       vjust=0),
        strip.text.x=element_text(face="bold"),
        strip.placement.y="outside",
        axis.text.y = element_text(hjust=-1),
        panel.spacing.x=unit(1,"lines"),
        panel.background = element_rect(fill = NA, color = "grey"),
        legend.position="none",
        panel.grid.major = element_line()
        )+
  labs(x="",y="")
fig1

ggsave(fig1,filename = "Figure1.PDF",dpi=300,width=10,height=6)


#############################################################
#### Figure 2: Estimate Results by Economic Liberalism  #####
#############################################################

CJa <- CJ[!is.na(CJ$ecoLR_bin),]

## Estimate Regression Model with interaction by Economic Left-Right
m2 <- setNames(lapply(att,function(x){
  svyglm(as.formula(paste("selected~",x,"*ecoLR_bin")),design=svydesign(data=CJa,id=~ResponseId))
}),att)

## Estimate Marginal Means
mm2 <-  data.frame(bind_rows(setNames(lapply(att,function(x){
  ggeffect(m2[[x]],terms=c(paste(x),"ecoLR_bin"))
}),att),.id="feature"))%>%
  mutate(group = str_replace(group,"ecoLR_bin",""))%>%
  rename(level=x,lower=conf.low,upper=conf.high,ecoLR=group)

## Estimate Average Marginal (Component) Effect by Sample
amce2 <- bind_rows(setNames(lapply(att,function(x){
  summary(margins::margins(m2[[x]],variables=paste(x),at=list(ecoLR_bin=levels(CJa$ecoLR_bin))))
}),att),.id="feature") %>% rename(level=factor,predicted=AME,ecoLR=ecoLR_bin) %>%
  mutate(.,level=str_replace(level,feature,""))

##Create Dummy values for reference categories
b <- bind_rows(lapply(att,function(x){
  expand.grid(colnames(CJa[paste(x)]), levels(CJa[[paste(x)]]),levels(CJa$ecoLR_bin))
}))
names(b) <- c("feature","level","ecoLR")
amce2 <- merge(b, amce2, by=c("feature","level","ecoLR"), all = T)
amce2[is.na(amce2)] <- 0

## Estimate Difference in Marginal Means between Samples
mmdiff2 <-  data.frame(bind_rows(setNames(lapply(att,function(x){
  summary(margins::margins(m2[[x]],variables=c("ecoLR_bin"),at=setNames(list(lvls[[x]]), paste(x))))
}),att),.id="feature")) %>%
  ungroup() %>%
  mutate_at(att,as.character) %>%
  mutate(level = as.factor(as.character(coalesce(!!! rlang::syms(att)))),
         factor = str_replace(factor,"ecoLR_bin","")) %>%
  select(.,-c(all_of(att))) %>% rename(predicted=AME,ecoLR=factor)


##Prepare for plotting
df2 <- bind_rows(mm2,amce2,mmdiff2,.id="ID")
levels(df2$level) <- lvl
df2$feature <- factor(df2$feature,levels=att,labels=ats)

df2$sig <- ifelse(
  (df2$ID == 1 & ((df2$lower < 0.5 & df2$upper < 0.5)|(df2$lower > 0.5 & df2$upper > 0.5)))|
    (df2$ID != 1 & ((df2$lower < 0 & df2$upper < 0)|(df2$lower > 0 & df2$upper > 0))),"not","yes")

df2$ID <- factor(df2$ID,levels=c(1,2,3),labels=c("Marginal Means","AMCE","Diff. Marginal Means"))

df2$level <- fct_relevel(df2$level,"Much more strict",after=20)
df2$level <- fct_relevel(df2$level,"Constant",after=20)

df2$ecoLR <- as.factor(as.character(df2$ecoLR))
df2$ecoLR <- fct_relevel(df2$ecoLR,"High",after=2)
##Plot Results by Sample

fig2 <- ggplot(df2, aes(x = level, y = predicted))+
  geom_hline(aes(yintercept=ifelse(ID=="Marginal Means",0.5,0)),size=1,color="grey")+
  geom_pointrange(aes(ymin=lower,
                      ymax=upper,
                      color=ecoLR,
                      shape=ecoLR,
                      alpha=sig),
                  position=position_dodge2(width=.75),
                  size=1.25,fatten=7)+
  facet_grid(feature~ID,
             scales="free",
             switch="y",
             space = "free")+
  scale_x_discrete(position = "top",
                   expand=c(-.65,0))+
  coord_flip()+
  theme_classic(base_family = "sans")+
  theme(text=element_text(size=16,family="Roboto Condensed"),
        strip.text.y.left=element_text(angle=0,hjust=0,
                                       vjust=0),
        strip.text.x=element_text(face="bold"),
        strip.placement.y="outside",
        axis.text.y = element_text(hjust=-1),
        panel.spacing.x=unit(1,"lines"),
        panel.background = element_rect(fill = NA, color = "grey"),
        legend.position="bottom",
        panel.grid.major = element_line(),
        legend.title=element_text(size=12))+
  labs(x="",y="")+
  scale_color_manual(values=c("#001524","#219ebc","#FF7D00"))+
  scale_alpha_manual(values=c(1,.25))+
  scale_shape_manual(values=c("\u25BC","\u25C6","\u25B2"))+
  guides(alpha="none",color=guide_legend(title="Economic Liberalism"),
         shape=guide_legend(title="Economic Liberalism"))
fig2

ggsave(fig2,filename="Figure2.PDF",width=14,height=8,dpi=500)


##############################################################
#### Figure 3: Estimate Results by Cultural Conservatism #####
##############################################################

CJa <- CJ[!is.na(CJ$culLR_bin),]

## Estimate Regression Model with interaction by Conservatism
m3 <- setNames(lapply(att,function(x){
  svyglm(as.formula(paste("selected~",x,"*culLR_bin")),design=svydesign(data=CJa,id=~ResponseId))
}),att)

## Estimate Marginal Means
mm3 <-  data.frame(bind_rows(setNames(lapply(att,function(x){
  ggeffect(m3[[x]],terms=c(paste(x),"culLR_bin"))
}),att),.id="feature"))%>%
  mutate(group = str_replace(group,"culLR_bin",""))%>%
  rename(level=x,lower=conf.low,upper=conf.high,culLR_bin=group)

## Estimate Average Marginal (Component) Effect by Conservatism
amce3 <- bind_rows(setNames(lapply(att,function(x){
  summary(margins::margins(m3[[x]],variables=paste(x),at=list(culLR_bin=levels(CJa$culLR_bin))))
}),att),.id="feature") %>% rename(level=factor,predicted=AME,culLR_bin=culLR_bin) %>%
  mutate(.,level=str_replace(level,feature,""))

##Create Dummy values for reference categories
b <- bind_rows(lapply(att,function(x){
  expand.grid(colnames(CJa[paste(x)]), levels(CJa[[paste(x)]]),levels(CJa$culLR_bin))
}))
names(b) <- c("feature","level","culLR_bin")
amce3 <- merge(b, amce3, by=c("feature","level","culLR_bin"), all = T)
amce3[is.na(amce3)] <- 0

## Estimate Difference in Marginal Means for Conservatism
mmdiff3 <-  data.frame(bind_rows(setNames(lapply(att,function(x){
  summary(margins::margins(m3[[x]],variables=c("culLR_bin"),at=setNames(list(lvls[[x]]), paste(x))))
}),att),.id="feature")) %>%
  ungroup() %>%
  mutate_at(att,as.character) %>%
  mutate(level = as.factor(as.character(coalesce(!!! rlang::syms(att)))),
         factor = str_replace(factor,"culLR_bin","")) %>%
  select(.,-c(all_of(att))) %>% rename(predicted=AME,culLR_bin=factor)


##Prepare for plotting
df3 <- bind_rows(mm3,amce3,mmdiff3,.id="ID")
levels(df3$level) <- lvl
df3$feature <- factor(df3$feature,levels=att,labels=ats)

df3$sig <- ifelse(
  (df3$ID == 1 & ((df3$lower < 0.5 & df3$upper < 0.5)|(df3$lower > 0.5 & df3$upper > 0.5)))|
    (df3$ID != 1 & ((df3$lower < 0 & df3$upper < 0)|(df3$lower > 0 & df3$upper > 0))),"not","yes")

df3$ID <- factor(df3$ID,levels=c(1,2,3),labels=c("Marginal Means","AMCE","Diff. Marginal Means"))

df3$level <- fct_relevel(df3$level,"Much more strict",after=20)
df3$level <- fct_relevel(df3$level,"Constant",after=20)

df3$culLR_bin <- as.factor(as.character(df3$culLR_bin))
df3$culLR_bin <- fct_relevel(df3$culLR_bin,"High",after=2)
##Plot Results by Conservatism

fig3 <- ggplot(df3, aes(x = level, y = predicted))+
  geom_hline(aes(yintercept=ifelse(ID=="Marginal Means",0.5,0)),size=1,color="grey")+
  geom_pointrange(aes(ymin=lower,
                      ymax=upper,
                      color=culLR_bin,
                      shape=culLR_bin,
                      alpha=sig),
                  position=position_dodge2(width=.75),
                  size=1.25,fatten=7)+
  facet_grid(feature~ID,
             scales="free",
             switch="y",
             space = "free")+
  scale_x_discrete(position = "top",
                   expand=c(-.65,0))+
  coord_flip()+
  theme_classic(base_family = "sans")+
  theme(text=element_text(size=16,family="Roboto Condensed"),
        strip.text.y.left=element_text(angle=0,hjust=0,
                                       vjust=0),
        strip.text.x=element_text(face="bold"),
        strip.placement.y="outside",
        axis.text.y = element_text(hjust=-1),
        panel.spacing.x=unit(1,"lines"),
        panel.background = element_rect(fill = NA, color = "grey"),
        legend.position="bottom",
        panel.grid.major = element_line(),
        legend.title=element_text(size=12))+
  labs(x="",y="")+
  scale_color_manual(values=c("#001524","#219ebc","#FF7D00"))+
  scale_alpha_manual(values=c(1,.25))+
  scale_shape_manual(values=c("\u25BC","\u25C6","\u25B2"))+
  guides(alpha="none",color=guide_legend(title="Cultural Conservatism"),
         shape=guide_legend(title="Cultural Conservatism"))
fig3

ggsave(fig3,filename="Figure3.PDF",width=14,height=7,dpi=500)


#####################################################################
#### Figure 4: Estimate Results by Attitudes towards Technology #####
#####################################################################

## Estimate Regression Model with interaction by Tech. Attitudes

m4 <- setNames(lapply(att,function(x){
  svyglm(as.formula(paste("selected~",x,"*techatt_bin")),design=svydesign(data=CJ,id=~ResponseId))
}),att)

## Estimate Marginal Means
mm4 <-  data.frame(bind_rows(setNames(lapply(att,function(x){
  ggeffect(m4[[x]],terms=c(paste(x),"techatt_bin"))
}),att),.id="feature"))%>%
  mutate(group = str_replace(group,"techatt_bin",""))%>%
  rename(level=x,lower=conf.low,upper=conf.high,tech=group)

## Estimate Average Marginal (Component) Effect by Sample
amce4 <- bind_rows(setNames(lapply(att,function(x){
  summary(margins::margins(m4[[x]],variables=paste(x),at=list(techatt_bin=levels(CJ$techatt_bin))))
}),att),.id="feature") %>% rename(level=factor,predicted=AME,tech=techatt_bin) %>%
  mutate(.,level=str_replace(level,feature,""))

##Create Dummy values for reference categories
b <- bind_rows(lapply(att,function(x){
  expand.grid(colnames(CJ[paste(x)]), levels(CJ[[paste(x)]]),levels(CJ$techatt_bin))
}))
names(b) <- c("feature","level","tech")
amce4 <- merge(b, amce4, by=c("feature","level","tech"), all = T)
amce4[is.na(amce4)] <- 0

## Estimate Difference in Marginal Means between Samples
mmdiff4 <-  data.frame(bind_rows(setNames(lapply(att,function(x){
  summary(margins::margins(m4[[x]],variables=c("techatt_bin"),at=setNames(list(lvls[[x]]), paste(x))))
}),att),.id="feature")) %>%
  ungroup() %>%
  mutate_at(att,as.character) %>%
  mutate(level = as.factor(as.character(coalesce(!!! rlang::syms(att)))),
         factor = str_replace(factor,"techatt_bin","")) %>%
  select(.,-c(all_of(att))) %>% rename(predicted=AME,tech=factor)


##Prepare for plotting
df4 <- bind_rows(mm4,amce4,mmdiff4,.id="ID")
levels(df4$level) <- lvl
df4$feature <- factor(df4$feature,levels=att,labels=ats)

df4$sig <- ifelse(
  (df4$ID == 1 & ((df4$lower < 0.5 & df4$upper < 0.5)|(df4$lower > 0.5 & df4$upper > 0.5)))|
    (df4$ID != 1 & ((df4$lower < 0 & df4$upper < 0)|(df4$lower > 0 & df4$upper > 0))),"not","yes")

df4$ID <- factor(df4$ID,levels=c(1,2,3),labels=c("Marginal Means","AMCE","Diff. Marginal Means"))

df4$level <- fct_relevel(df4$level,"Much more strict",after=20)
df4$level <- fct_relevel(df4$level,"Constant",after=20)

df4$tech <- as.factor(as.character(df4$tech))
df4$tech <- fct_relevel(df4$tech,"High",after=2)
##Plot Results by Sample

fig4 <- ggplot(df4, aes(x = level, y = predicted))+
  geom_hline(aes(yintercept=ifelse(ID=="Marginal Means",0.5,0)),linewidth=1,color="grey")+
  geom_pointrange(aes(ymin=lower,
                      ymax=upper,
                      color=tech,
                      shape=tech,
                      alpha=sig),
                  position=position_dodge2(width=.75),
                  size=1.25,fatten=7)+
  facet_grid(feature~ID,
             scales="free",
             switch="y",
             space = "free")+
  scale_x_discrete(position = "top",
                   expand=c(-.65,0))+
  scale_y_continuous(breaks=seq(-1,1,by=.05))+
  coord_flip()+
  theme_classic(base_family = "sans")+
  theme(text=element_text(size=16,family="Roboto Condensed"),
        strip.text.y.left=element_text(angle=0,hjust=0,
                                       vjust=0),
        strip.text.x=element_text(face="bold"),
        strip.placement.y="outside",
        axis.text.y = element_text(hjust=-1),
        panel.spacing.x=unit(1,"lines"),
        panel.background = element_rect(fill = NA, color = "grey"),
        legend.position="bottom",
        panel.grid.major = element_line(),
        legend.title=element_text(size=12))+
  labs(x="",y="")+
  scale_color_manual(values=c("#001524","#219ebc","#FF7D00"))+
  scale_alpha_manual(values=c(1,.25))+
  scale_shape_manual(values=c("\u25BC","\u25C6","\u25B2"))+
  guides(alpha="none",color=guide_legend(title="Technology Optimism"),
         shape=guide_legend(title="Technology Optimism"))
fig4

ggsave(fig4,filename="Figure4.PDF",width=14,height=8,dpi=500)

#####################################
#### Figure 5:  Results by Sample #####
#####################################

## Estimate Regression Model with interaction by Sample
CJ$sample <- as.factor(CJ$sample)

m5 <- setNames(lapply(att,function(x){
  svyglm(as.formula(paste("selected~",x,"*sample")),design=svydesign(data=CJ,id=~ResponseId))
}),att)

## Estimate Marginal Means by Sample
mm5 <-  data.frame(bind_rows(setNames(lapply(att,function(x){
  as.data.frame(ggeffect(m5[[x]],terms=c(paste(x),"sample")))
}),att),.id="feature"))%>%
  mutate(group = str_replace(group,"sample",""))%>%
  rename(level=x,lower=conf.low,upper=conf.high,sample=group)

## Estimate Average Marginal (Component) Effect by Sample
amce5 <- bind_rows(setNames(lapply(att,function(x){
  as.data.frame(summary(margins::margins(m5[[x]],variables=paste(x),at=list(sample=levels(CJ$sample)))))
}),att),.id="feature") %>% rename(level=factor,predicted=AME) %>%
  mutate(.,level=str_replace(level,feature,""))

##Create Dummy values for reference categories
b <- bind_rows(lapply(att,function(x){
  expand.grid(colnames(CJ[paste(x)]), levels(CJ[[paste(x)]]),levels(CJ$sample))
}))
names(b) <- c("feature","level","sample")
amce5 <- merge(b, amce5, by=c("feature","level","sample"), all = T)
amce5[is.na(amce5)] <- 0

## Estimate Difference in Marginal Means between Samples
mmdiff5 <-  data.frame(bind_rows(setNames(lapply(att,function(x){
  as.data.frame(summary(margins::margins(m5[[x]],variables=c("sample"),at=setNames(list(lvls[[x]]), paste(x)))))
}),att),.id="feature")) %>%
  ungroup() %>%
  mutate_at(att,as.character) %>%
  mutate(level = as.factor(as.character(coalesce(!!! rlang::syms(att)))),
         factor = str_replace(factor,"sample","")) %>%
  select(.,-c(all_of(att))) %>% rename(predicted=AME,sample=factor)

### Change reference for Diff. in MM

CJ$sample <- fct_rev(CJ$sample)

m5a <- setNames(lapply(att,function(x){
  svyglm(as.formula(paste("selected~",x,"*sample")),design=svydesign(data=CJ,id=~ResponseId))
}),att)

mmdiff5a <-  data.frame(bind_rows(setNames(lapply(att,function(x){
  summary(margins::margins(m5a[[x]],variables=c("sample"),at=setNames(list(lvls[[x]]), paste(x))))
}),att),.id="feature")) %>%
  ungroup() %>%
  mutate_at(att,as.character) %>%
  mutate(level = as.factor(as.character(coalesce(!!! rlang::syms(att)))),
         factor = str_replace(factor,"sample","")) %>%
  select(.,-c(all_of(att))) %>% rename(predicted=AME,sample=factor)


##Prepare for plotting
df5 <- bind_rows(mm5,amce5,mmdiff5,mmdiff5a[mmdiff5a$sample=="JP",],.id="ID")

levels(df5$level) <- lvl
df5$feature <- factor(df5$feature,levels=att,labels=ats)

df5$sig <- ifelse(
  (df5$ID == 1 & ((df5$lower < 0.5 & df5$upper < 0.5)|(df5$lower > 0.5 & df5$upper > 0.5)))|
    (df5$ID != 1 & ((df5$lower < 0 & df5$upper < 0)|(df5$lower > 0 & df5$upper > 0))),"not","yes")

df5$reference <- case_when(df5$ID==3 ~ "vs. DE",
                           df5$ID==4 ~ "vs. US")

df5$ID <- factor(df5$ID,levels=c(1,2,3,4),labels=c("Marginal Means","AMCE","Diff. Marginal Means\nref:DE (dash US)","Diff. Marginal Means\nref:DE (dash US)"))


df5$level <- fct_relevel(df5$level,"Much more strict",after=20)
df5$level <- fct_relevel(df5$level,"Constant",after=20)

##Plot Results by Sample

fig5 <- ggplot(df5, aes(x = level, y = predicted))+
  geom_hline(aes(yintercept=ifelse(ID=="Marginal Means",0.5,0)),size=1,color="grey")+
  geom_pointrange(aes(ymin=lower,
                      ymax=upper,
                      color=sample,
                      shape=sample,
                      alpha=sig,
                      linetype=ifelse(!is.na(reference),reference,"solid"))
                  ,position=position_dodge2(width=.75),
                  size=1.25,fatten=3)+
  facet_grid(feature~ID,
             scales="free",
             switch="y",
             space = "free")+
  scale_x_discrete(position = "top",
                   expand=c(-.65,0))+
  scale_y_continuous(breaks=seq(-1,1,by=.05))+
  coord_flip()+
  theme_classic(base_family = "sans")+
  theme(text=element_text(size=16,family="Roboto Condensed"),
        strip.text.y.left=element_text(angle=0,hjust=0,
                                       vjust=0),
        strip.text.x=element_text(face="bold"),
        strip.placement.y="outside",
        axis.text.y = element_text(hjust=-1),
        panel.spacing.x=unit(1,"lines"),
        panel.background = element_rect(fill = NA, color = "grey"),
        legend.position="bottom",
        panel.grid.major = element_line(),
        legend.title=element_blank())+
  labs(x="",y="")+
  scale_color_manual(values=c("#ffa700","#d62d20","navy","#ffa700","#d62d20","navy"))+
  scale_alpha_manual(values=c(1,.5))+
  scale_shape_manual(values=c(15,16,17))+
  scale_linetype_manual(values=c("solid","solid","dashed"))+
  guides(alpha="none",linetype="none")
fig5

ggsave(fig5,filename="Figure5.PDF",width=12,height=8,dpi=500)


###################
#### Appendix #####
###################

################################################################
#### Appendix Figure 1: Distribution of Attitudinal Measures ###
################################################################

dfa1 <- rbind(
  select(CJ,c(ecoLR,sample,ecoLR_bin)) %>%
    rename("x"=1,"sample"=2,"bin"=3) %>%
    mutate(var = "Economic Liberalism"),
  select(CJ,c(culLR,sample,culLR_bin)) %>%
    rename("x"=1,"sample"=2,"bin"=3) %>%
    mutate(var = "Cultural Conservatism"),
  select(CJ,c(techatt,sample,techatt_bin)) %>%
    rename("x"=1,"sample"=2,"bin"=3) %>%
    mutate(var = "Technology Optimism")
  ) %>% mutate(var = fct_relevel(as.factor(var),"Economic Liberalism",after=0),
               bin = as.factor(bin))

pa1 <- ggplot(dfa1,aes(x=x,y=sample,fill=sample,linetype=bin))+
  geom_density_ridges2(scale=2,bandwidth=0.5,alpha=0.9,size=.7)+
  theme_classic()+
  scale_fill_manual(values=c("#F14A16","#370665","#35589A"))+
  labs(y="",x="")+
  guides(fill="none")+
  facet_wrap(~var,scales="free_x",ncol=3)+
  theme(text=element_text(size=16),
        legend.title=element_blank(),
        legend.position="bottom")+
  scale_x_continuous(breaks=seq(1,10,by=1))
pa1
ggsave(pa1,filename="Appendix_Figure1.PDF",dpi=500,width=8,height=5)


################################################
##### Appendix Figure3: Carryover Effects   ####
################################################

CJa <- CJ
CJa$task <- as.factor(as.character(CJa$task))

## Estimate Regression Model with interaction by task
ma3 <- setNames(lapply(att,function(x){
  svyglm(as.formula(paste("selected~",x,"*task")),design=svydesign(data=CJa,id=~ResponseId))
}),att)
## Estimate Marginal Means
mma3 <-  data.frame(bind_rows(setNames(lapply(att,function(x){
  ggeffect(ma3[[x]],terms=c(paste(x),"task"))
}),att),.id="feature"))%>%
  mutate(group = str_replace(group,"task",""))%>%
  rename(level=x,lower=conf.low,upper=conf.high,task=group)
## Estimate Average Marginal (Component) Effect by task
amcea3 <- bind_rows(setNames(lapply(att,function(x){
  summary(margins::margins(ma3[[x]],variables=paste(x),at=list(task=levels(CJa$task))))
}),att),.id="feature") %>% rename(level=factor,predicted=AME,task=task) %>%
  mutate(.,level=str_replace(level,feature,""))
##Create Dummy values for reference categories
b <- bind_rows(lapply(att,function(x){
  expand.grid(colnames(CJa[paste(x)]), levels(CJa[[paste(x)]]),levels(CJa$task))
}))
names(b) <- c("feature","level","task")
amcea3 <- merge(b, amcea3, by=c("feature","level","task"), all = T)
amcea3[is.na(amcea3)] <- 0
## Estimate Difference in Marginal Means between task
mmdiffa3 <-  data.frame(bind_rows(setNames(lapply(att,function(x){
  summary(margins::margins(ma3[[x]],variables=c("task"),at=setNames(list(lvls[[x]]), paste(x))))
}),att),.id="feature")) %>%
  ungroup() %>%
  mutate_at(att,as.character) %>%
  mutate(level = as.factor(as.character(coalesce(!!! rlang::syms(att)))),
         factor = str_replace(factor,"task","")) %>%
  select(.,-c(all_of(att))) %>% rename(predicted=AME,task=factor)
##Prepare for plotting
dfa3 <- bind_rows(mma3,amcea3,mmdiffa3,.id="ID")
levels(dfa3$level) <- lvl
dfa3$feature <- factor(dfa3$feature,levels=att,labels=ats)

dfa3$sig <- ifelse(
  (dfa3$ID == 1 & ((dfa3$lower < 0.5 & dfa3$upper < 0.5)|(dfa3$lower > 0.5 & dfa3$upper > 0.5)))|
    (dfa3$ID != 1 & ((dfa3$lower < 0 & dfa3$upper < 0)|(dfa3$lower > 0 & dfa3$upper > 0))),"not","yes")

dfa3$ID <- factor(dfa3$ID,levels=c(1,2,3),labels=c("Marginal Means","AMCE","Diff. Marginal Means"))

dfa3$level <- fct_relevel(dfa3$level,"Much more strict",after=20)
dfa3$level <- fct_relevel(dfa3$level,"Constant",after=20)

dfa3$task <- as.factor(as.character(dfa3$task))
dfa3$task <- fct_relevel(dfa3$task,"High",after=2)

##Plot Results by Task

figa3 <- ggplot(dfa3, aes(x = level, y = predicted))+
  geom_hline(aes(yintercept=ifelse(ID=="Marginal Means",0.5,0)),size=1,color="grey")+
  geom_pointrange(aes(ymin=lower,
                      ymax=upper,
                      color=task,
                      shape=task,
                      alpha=sig),
                  position=position_dodge2(width=.75),
                  size=1.25,fatten=1)+
  facet_grid(feature~ID,
             scales="free",
             switch="y",
             space = "free")+
  scale_x_discrete(position = "top",
                   expand=c(-.65,0))+
  coord_flip()+
  theme_classic(base_family = "sans")+
  theme(text=element_text(size=16,family="Roboto Condensed"),
        strip.text.y.left=element_text(angle=0,hjust=0,
                                       vjust=0),
        strip.text.x=element_text(face="bold"),
        strip.placement.y="outside",
        axis.text.y = element_text(hjust=-1),
        panel.spacing.x=unit(1,"lines"),
        panel.background = element_rect(fill = NA, color = "grey"),
        legend.position="bottom",
        panel.grid.major = element_line(),
        legend.title=element_text(size=12))+
  labs(x="",y="")+
  scale_color_manual(values=c("#005f60","#008083","#249ea0","#faab36","#f78104","#fd5901"))+
  scale_alpha_manual(values=c(1,.25))+
  scale_shape_manual(values=c(rep(15,9)))+
  guides(alpha="none",color=guide_legend(title="task"),
         shape=guide_legend(title="task"))
figa3

ggsave(figa3,filename="Appendix_Figure3.PDF",width=12,height=7,dpi=500)

#####################################################
##### Appendix Figure 4: Profile Order Effects   ####
#####################################################

CJa <- CJ
CJa$profile <- as.factor(as.character(CJa$profile))

## Estimate Regression Model with interaction by profile
ma4 <- setNames(lapply(att,function(x){
  svyglm(as.formula(paste("selected~",x,"*profile")),design=svydesign(data=CJa,id=~ResponseId))
}),att)
## Estimate Marginal Means
mma4 <-  data.frame(bind_rows(setNames(lapply(att,function(x){
  ggeffect(ma4[[x]],terms=c(paste(x),"profile"))
}),att),.id="feature"))%>%
  mutate(group = str_replace(group,"profile",""))%>%
  rename(level=x,lower=conf.low,upper=conf.high,profile=group)
## Estimate Average Marginal (Component) Effect by profile
amcea4 <- bind_rows(setNames(lapply(att,function(x){
  summary(margins::margins(ma4[[x]],variables=paste(x),at=list(profile=levels(CJa$profile))))
}),att),.id="feature") %>% rename(level=factor,predicted=AME,profile=profile) %>%
  mutate(.,level=str_replace(level,feature,""))
##Create Dummy values for reference categories
b <- bind_rows(lapply(att,function(x){
  expand.grid(colnames(CJa[paste(x)]), levels(CJa[[paste(x)]]),levels(CJa$profile))
}))
names(b) <- c("feature","level","profile")
amcea4 <- merge(b, amcea4, by=c("feature","level","profile"), all = T)
amcea4[is.na(amcea4)] <- 0
## Estimate Difference in Marginal Means between profile
mmdiffa4 <-  data.frame(bind_rows(setNames(lapply(att,function(x){
  summary(margins::margins(ma4[[x]],variables=c("profile"),at=setNames(list(lvls[[x]]), paste(x))))
}),att),.id="feature")) %>%
  ungroup() %>%
  mutate_at(att,as.character) %>%
  mutate(level = as.factor(as.character(coalesce(!!! rlang::syms(att)))),
         factor = str_replace(factor,"profile","")) %>%
  select(.,-c(all_of(att))) %>% rename(predicted=AME,profile=factor)
##Prepare for plotting
dfa4 <- bind_rows(mma4,amcea4,mmdiffa4,.id="ID")
levels(dfa4$level) <- lvl
dfa4$feature <- factor(dfa4$feature,levels=att,labels=ats)

dfa4$sig <- ifelse(
  (dfa4$ID == 1 & ((dfa4$lower < 0.5 & dfa4$upper < 0.5)|(dfa4$lower > 0.5 & dfa4$upper > 0.5)))|
    (dfa4$ID != 1 & ((dfa4$lower < 0 & dfa4$upper < 0)|(dfa4$lower > 0 & dfa4$upper > 0))),"not","yes")

dfa4$ID <- factor(dfa4$ID,levels=c(1,2,3),labels=c("Marginal Means","AMCE","Diff. Marginal Means"))

dfa4$level <- fct_relevel(dfa4$level,"Much more strict",after=20)
dfa4$level <- fct_relevel(dfa4$level,"Constant",after=20)

dfa4$profile <- as.factor(as.character(dfa4$profile))
dfa4$profile <- fct_relevel(dfa4$profile,"High",after=2)

##Plot Results by Profile Order

figa4 <- ggplot(dfa4, aes(x = level, y = predicted))+
  geom_hline(aes(yintercept=ifelse(ID=="Marginal Means",0.5,0)),size=1,color="grey")+
  geom_pointrange(aes(ymin=lower,
                      ymax=upper,
                      color=profile,
                      shape=profile,
                      alpha=sig),
                  position=position_dodge2(width=.75),
                  size=1.25,fatten=1)+
  facet_grid(feature~ID,
             scales="free",
             switch="y",
             space = "free")+
  scale_x_discrete(position = "top",
                   expand=c(-.65,0))+
  coord_flip()+
  theme_classic(base_family = "sans")+
  theme(text=element_text(size=16,family="Roboto Condensed"),
        strip.text.y.left=element_text(angle=0,hjust=0,
                                       vjust=0),
        strip.text.x=element_text(face="bold"),
        strip.placement.y="outside",
        axis.text.y = element_text(hjust=-1),
        panel.spacing.x=unit(1,"lines"),
        panel.background = element_rect(fill = NA, color = "grey"),
        legend.position="bottom",
        panel.grid.major = element_line(),
        legend.title=element_text(size=12))+
  labs(x="",y="")+
  scale_color_manual(values=c("#005f60","#fd5901"))+
  scale_alpha_manual(values=c(1,.25))+
  scale_shape_manual(values=c(rep(15,9)))+
  guides(alpha="none",color=guide_legend(title="profile"),
         shape=guide_legend(title="profile"))
figa4

ggsave(figa4,filename="Appendix_Figure4.PDF",width=12,height=7,dpi=500)


#############################################
##### Appendix Figure 5: Balancing Tests ####
#############################################

CJa <- CJ
CJa$female <- ifelse(CJa$gender == "Male",0,1)
CJa$age <- as.numeric(as.character(CJa$age))
CJa$sample <- as.numeric(as.factor(CJa$sample))

ma5a <- setNames(lapply(att,function(x){
  svyglm(as.formula(paste("ecoLR~",x)),design=svydesign(data=CJa,id=~ResponseId))
}),att)
ma5b <- setNames(lapply(att,function(x){
  svyglm(as.formula(paste("culLR~",x)),design=svydesign(data=CJa,id=~ResponseId))
}),att)
ma5c <- setNames(lapply(att,function(x){
  svyglm(as.formula(paste("techatt~",x)),design=svydesign(data=CJa,id=~ResponseId))
}),att)
ma5d <- setNames(lapply(att,function(x){
  svyglm(as.formula(paste("sample~",x)),design=svydesign(data=CJa,id=~ResponseId))
}),att)
ma5e <- setNames(lapply(att,function(x){
  svyglm(as.formula(paste("female~",x)),design=svydesign(data=CJa,id=~ResponseId))
}),att)
ma5f <- setNames(lapply(att,function(x){
  svyglm(as.formula(paste("age~",x)),design=svydesign(data=CJa,id=~ResponseId))
}),att)

dfa5 <- bind_rows(data.frame(bind_rows(setNames(lapply(att,function(x){
  ggeffect(ma5a[[x]],paste(x))
}),att),.id="feature")) %>% mutate(.,vline=mean(CJa$ecoLR[!is.na(CJa$ecoLR)])),
data.frame(bind_rows(setNames(lapply(att,function(x){
  ggeffect(ma5b[[x]],paste(x))
}),att),.id="feature"))%>% mutate(.,vline=mean(CJa$culLR[!is.na(CJa$culLR)])),
data.frame(bind_rows(setNames(lapply(att,function(x){
  ggeffect(ma5c[[x]],paste(x))
}),att),.id="feature"))%>% mutate(.,vline=mean(CJa$techatt[!is.na(CJa$techatt)])),
data.frame(bind_rows(setNames(lapply(att,function(x){
  ggeffect(ma5d[[x]],paste(x))
}),att),.id="feature"))%>% mutate(.,vline=mean(CJa$sample[!is.na(CJa$sample)])),
data.frame(bind_rows(setNames(lapply(att,function(x){
  ggeffect(ma5e[[x]],paste(x))
}),att),.id="feature"))%>% mutate(.,vline=mean(CJa$female[!is.na(CJa$female)])),
data.frame(bind_rows(setNames(lapply(att,function(x){
  ggeffect(ma5f[[x]],paste(x))
}),att),.id="feature"))%>% mutate(.,vline=mean(CJa$age[!is.na(CJa$age)])),.id="ID")

levels(dfa5$x) <- lvl
dfa5$feature <- factor(dfa5$feature,levels=att,labels=ats)

dfa5$sig <- ifelse((dfa5$conf.low < dfa5$vline & dfa5$conf.high < dfa5$vline)|(dfa5$conf.low > dfa5$vline & dfa5$conf.high > dfa5$vline),"not","yes")
dfa5$ID <- factor(dfa5$ID,levels=c(1,2,3,4,5,6),labels=c("Economic Liberalism",
                                                         "Cultural Conservatism",
                                                         "Tech. Optimism",
                                                         "Sample",
                                                         "Female vs. Male",
                                                         "Age"))

dfa5$x <- fct_relevel(dfa5$x,"Much more strict",after=20)
dfa5$x <- fct_relevel(dfa5$x,"Constant",after=20)

figa5 <- ggplot(dfa5, aes(x = x, y = predicted))+
  geom_hline(aes(yintercept=vline),size=1,color="grey")+
  geom_pointrange(aes(ymin=conf.low,
                      ymax=conf.high,
                      alpha=sig),
                  position=position_dodge2(width=.75),
                  size=1.25,fatten=1,color="blue")+
  facet_grid(feature~ID,
             scales="free",
             switch="y")+
  scale_x_discrete(position = "top")+
  coord_flip()+
  theme_classic(base_family = "sans")+
  theme(text=element_text(size=16,family="Roboto Condensed"),
        strip.text.y.left=element_text(angle=0,hjust=0,
                                       vjust=0),
        strip.text.x=element_text(face="bold"),
        strip.placement.y="outside",
        axis.text.y = element_text(hjust=-1),
        panel.spacing.x=unit(1,"lines"),
        panel.background = element_rect(fill = NA, color = "grey"),
        legend.position="bottom",
        panel.grid.major = element_line(),
        legend.title=element_text(size=12))+
  labs(x="",y="")+
  scale_color_manual(values=c("#005f60","#fd5901"))+
  scale_alpha_manual(values=c(1,.25))+
  scale_shape_manual(values=c(rep(15,9)))+
  guides(alpha="none",color=guide_legend(title="profile"),
         shape=guide_legend(title="profile"))
figa5

ggsave(figa5,filename="Appendix_Figure5.PDF",width=15,height=8,dpi=500)

########################################################################
##### Appendix Figure 6a: Results by Voting Intention (United States) ##
########################################################################

CJ$vote <- as.factor(as.character(CJ$vote))
levels(CJ$vote)

## US
CJa <- CJ[CJ$sample=="US",]
CJa <- CJa[CJa$vote == "Republicans" | CJa$vote == "Democrats",]
CJa$vote <- as.factor(as.character(CJa$vote))

## Estimate Regression Model with interaction by vote
ma6a <- setNames(lapply(att,function(x){
  svyglm(as.formula(paste("selected~",x,"*vote")),design=svydesign(data=CJa,id=~ResponseId))
}),att)
## Estimate Marginal Means
mma6a <-  data.frame(bind_rows(setNames(lapply(att,function(x){
  ggeffect(ma6a[[x]],terms=c(paste(x),"vote"))
}),att),.id="feature"))%>%
  mutate(group = str_replace(group,"vote",""))%>%
  rename(level=x,lower=conf.low,upper=conf.high,vote=group)
## Estimate Average Marginal (Component) Effect by vote
amcea6a <- bind_rows(setNames(lapply(att,function(x){
  summary(margins::margins(ma6a[[x]],variables=paste(x),at=list(vote=levels(CJa$vote))))
}),att),.id="feature") %>% rename(level=factor,predicted=AME,vote=vote) %>%
  mutate(.,level=str_replace(level,feature,""))
##Create Dummy values for reference categories
b <- bind_rows(lapply(att,function(x){
  expand.grid(colnames(CJa[paste(x)]), levels(CJa[[paste(x)]]),levels(CJa$vote))
}))
names(b) <- c("feature","level","vote")
amcea6a <- merge(b, amcea6a, by=c("feature","level","vote"), all = T)
amcea6a[is.na(amcea6a)] <- 0
## Estimate Difference in Marginal Means between vote
mmdiffa6a <-  data.frame(bind_rows(setNames(lapply(att,function(x){
  summary(margins::margins(ma6a[[x]],variables=c("vote"),at=setNames(list(lvls[[x]]), paste(x))))
}),att),.id="feature")) %>%
  ungroup() %>%
  mutate_at(att,as.character) %>%
  mutate(level = as.factor(as.character(coalesce(!!! rlang::syms(att)))),
         factor = str_replace(factor,"vote","")) %>%
  select(.,-c(all_of(att))) %>% rename(predicted=AME,vote=factor)
##Prepare for plotting
dfa6a <- bind_rows(mma6a,amcea6a,mmdiffa6a,.id="ID")
levels(dfa6a$level) <- lvl
dfa6a$feature <- factor(dfa6a$feature,levels=att,labels=ats)

dfa6a$sig <- ifelse(
  (dfa6a$ID == 1 & ((dfa6a$lower < 0.5 & dfa6a$upper < 0.5)|(dfa6a$lower > 0.5 & dfa6a$upper > 0.5)))|
    (dfa6a$ID != 1 & ((dfa6a$lower < 0 & dfa6a$upper < 0)|(dfa6a$lower > 0 & dfa6a$upper > 0))),"not","yes")

dfa6a$ID <- factor(dfa6a$ID,levels=c(1,2,3),labels=c("Marginal Means","AMCE","Diff. Marginal Means"))

dfa6a$level <- fct_relevel(dfa6a$level,"Much more strict",after=20)
dfa6a$level <- fct_relevel(dfa6a$level,"Constant",after=20)

dfa6a$vote <- as.factor(as.character(dfa6a$vote))

##Plot Results by Voting Intention

figa6a <- ggplot(dfa6a, aes(x = level, y = predicted))+
  geom_hline(aes(yintercept=ifelse(ID=="Marginal Means",0.5,0)),size=1,color="grey")+
  geom_pointrange(aes(ymin=lower,
                      ymax=upper,
                      color=vote,
                      shape=vote,
                      alpha=sig),
                  position=position_dodge2(width=.75),
                  size=1.25,fatten=1)+
  facet_grid(feature~ID,
             scales="free",
             switch="y",
             space = "free")+
  scale_x_discrete(position = "top",
                   expand=c(-.65,0))+
  coord_flip()+
  theme_classic(base_family = "sans")+
  theme(text=element_text(size=16,family="Roboto Condensed"),
        strip.text.y.left=element_text(angle=0,hjust=0,
                                       vjust=0),
        strip.text.x=element_text(face="bold"),
        strip.placement.y="outside",
        axis.text.y = element_text(hjust=-1),
        panel.spacing.x=unit(1,"lines"),
        panel.background = element_rect(fill = NA, color = "grey"),
        legend.position="bottom",
        panel.grid.major = element_line(),
        legend.title=element_text(size=12))+
  labs(x="",y="")+
  scale_color_manual(values=c("blue","red"))+
  scale_alpha_manual(values=c(1,.25))+
  scale_shape_manual(values=c(rep(15,9)))+
  guides(alpha="none",color=guide_legend(title="vote"),
         shape=guide_legend(title="vote"))
figa6a

ggsave(figa6a,filename="Appendix_Figure6a.PDF",width=12,height=7,dpi=500)

########################################################################
##### Appendix Figure 6b: Results by Voting Intention (Germany) ########
########################################################################

CJa <- CJ[CJ$sample=="DE",]
CJa <- CJa[CJa$vote != "Would not vote" & CJa$vote != "Other",]
CJa$vote <- as.factor(as.character(CJa$vote))
CJa$vote <- fct_rev(CJa$vote)
levels(CJa$vote)
CJa$vote <- fct_relevel(CJa$vote,"SPD",after=2)

## Estimate Regression Model with interaction by vote
ma6b <- setNames(lapply(att,function(x){
  svyglm(as.formula(paste("selected~",x,"*vote")),design=svydesign(data=CJa,id=~ResponseId))
}),att)
## Estimate Marginal Means
mma6b <-  data.frame(bind_rows(setNames(lapply(att,function(x){
  ggeffect(ma6b[[x]],terms=c(paste(x),"vote"))
}),att),.id="feature"))%>%
  mutate(group = str_replace(group,"vote",""))%>%
  rename(level=x,lower=conf.low,upper=conf.high,vote=group)
## Estimate Average Marginal (Component) Effect by vote
amcea6b <- bind_rows(setNames(lapply(att,function(x){
  summary(margins::margins(ma6b[[x]],variables=paste(x),at=list(vote=levels(CJa$vote))))
}),att),.id="feature") %>% rename(level=factor,predicted=AME,vote=vote) %>%
  mutate(.,level=str_replace(level,feature,""))
##Create Dummy values for reference categories
b <- bind_rows(lapply(att,function(x){
  expand.grid(colnames(CJa[paste(x)]), levels(CJa[[paste(x)]]),levels(CJa$vote))
}))
names(b) <- c("feature","level","vote")
amcea6b <- merge(b, amcea6b, by=c("feature","level","vote"), all = T)
amcea6b[is.na(amcea6b)] <- 0
## Estimate Difference in Marginal Means between vote
mmdiffa6b <-  data.frame(bind_rows(setNames(lapply(att,function(x){
  summary(margins::margins(ma6b[[x]],variables=c("vote"),at=setNames(list(lvls[[x]]), paste(x))))
}),att),.id="feature")) %>%
  ungroup() %>%
  mutate_at(att,as.character) %>%
  mutate(level = as.factor(as.character(coalesce(!!! rlang::syms(att)))),
         factor = str_replace(factor,"vote","")) %>%
  select(.,-c(all_of(att))) %>% rename(predicted=AME,vote=factor)
##Prepare for plotting
dfa6b <- bind_rows(mma6b,amcea6b,mmdiffa6b,.id="ID")
levels(dfa6b$level) <- lvl
dfa6b$feature <- factor(dfa6b$feature,levels=att,labels=ats)

dfa6b$sig <- ifelse(
  (dfa6b$ID == 1 & ((dfa6b$lower < 0.5 & dfa6b$upper < 0.5)|(dfa6b$lower > 0.5 & dfa6b$upper > 0.5)))|
    (dfa6b$ID != 1 & ((dfa6b$lower < 0 & dfa6b$upper < 0)|(dfa6b$lower > 0 & dfa6b$upper > 0))),"not","yes")

dfa6b$ID <- factor(dfa6b$ID,levels=c(1,2,3),labels=c("Marginal Means","AMCE","Diff. Marginal Means"))

dfa6b$level <- fct_relevel(dfa6b$level,"Much more strict",after=20)
dfa6b$level <- fct_relevel(dfa6b$level,"Constant",after=20)

dfa6b$vote <- as.factor(as.character(dfa6b$vote))
levels(dfa6b$vote)
dfa6b$vote <- fct_rev(dfa6b$vote)
dfa6b$vote <- fct_relevel(dfa6b$vote,"SPD",after=2)
dfa6b$vote <- fct_relevel(dfa6b$vote,"FDP",after=4)

##Plot Results by Voting Intention

figa6b <- ggplot(dfa6b, aes(x = level, y = predicted))+
  geom_hline(aes(yintercept=ifelse(ID=="Marginal Means",0.5,0)),size=1,color="grey")+
  geom_pointrange(aes(ymin=lower,
                      ymax=upper,
                      color=vote,
                      shape=vote,
                      alpha=sig),
                  position=position_dodge2(width=.75),
                  size=1.25,fatten=1)+
  facet_grid(feature~ID,
             scales="free",
             switch="y",
             space = "free")+
  scale_x_discrete(position = "top",
                   expand=c(-.65,0))+
  coord_flip()+
  theme_classic(base_family = "sans")+
  theme(text=element_text(size=16,family="Roboto Condensed"),
        strip.text.y.left=element_text(angle=0,hjust=0,
                                       vjust=0),
        strip.text.x=element_text(face="bold"),
        strip.placement.y="outside",
        axis.text.y = element_text(hjust=-1),
        panel.spacing.x=unit(1,"lines"),
        panel.background = element_rect(fill = NA, color = "grey"),
        legend.position="bottom",
        panel.grid.major = element_line(),
        legend.title=element_text(size=12))+
  labs(x="",y="")+
  scale_color_manual(values=c("#A6006B","#1AA037","#E3000F","#000000","#FFEF00","#0489DB"))+
  scale_alpha_manual(values=c(1,.25))+
  scale_shape_manual(values=c(rep(15,9)))+
  guides(alpha="none",color=guide_legend(title="vote"),
         shape=guide_legend(title="vote"))
figa6b

ggsave(figa6b,filename="Appendix_Figure6b.PDF",width=12,height=7,dpi=500)


##############################################################
##### Appendix Figure 7: (Tech7) Technology solves issues ####
##############################################################

CJa <- CJ[!is.na(CJ$tech7),]
CJa$tech7 <- factor(CJa$tech7,levels=c(1,2,3,4,5),labels=c("Low","Low","Low","High","High"))

## Estimate Regression Model with interaction by tech7
ma7 <- setNames(lapply(att,function(x){
  svyglm(as.formula(paste("selected~",x,"*tech7")),design=svydesign(data=CJa,id=~ResponseId))
}),att)
## Estimate Marginal Means
mma7 <-  data.frame(bind_rows(setNames(lapply(att,function(x){
  ggeffect(ma7[[x]],terms=c(paste(x),"tech7"))
}),att),.id="feature"))%>%
  mutate(group = str_replace(group,"tech7",""))%>%
  rename(level=x,lower=conf.low,upper=conf.high,tech7=group)
## Estimate Average Marginal (Component) Effect by tech7
amcea7 <- bind_rows(setNames(lapply(att,function(x){
  summary(margins::margins(ma7[[x]],variables=paste(x),at=list(tech7=levels(CJa$tech7))))
}),att),.id="feature") %>% rename(level=factor,predicted=AME,tech7=tech7) %>%
  mutate(.,level=str_replace(level,feature,""))
##Create Dummy values for reference categories
b <- bind_rows(lapply(att,function(x){
  expand.grid(colnames(CJa[paste(x)]), levels(CJa[[paste(x)]]),levels(CJa$tech7))
}))
names(b) <- c("feature","level","tech7")
amcea7 <- merge(b, amcea7, by=c("feature","level","tech7"), all = T)
amcea7[is.na(amcea7)] <- 0
## Estimate Difference in Marginal Means between tech7
mmdiffa7 <-  data.frame(bind_rows(setNames(lapply(att,function(x){
  summary(margins::margins(ma7[[x]],variables=c("tech7"),at=setNames(list(lvls[[x]]), paste(x))))
}),att),.id="feature")) %>%
  ungroup() %>%
  mutate_at(att,as.character) %>%
  mutate(level = as.factor(as.character(coalesce(!!! rlang::syms(att)))),
         factor = str_replace(factor,"tech7","")) %>%
  select(.,-c(all_of(att))) %>% rename(predicted=AME,tech7=factor)
##Prepare for plotting
dfa7 <- bind_rows(mma7,amcea7,mmdiffa7,.id="ID")
levels(dfa7$level) <- lvl
dfa7$feature <- factor(dfa7$feature,levels=att,labels=ats)

dfa7$sig <- ifelse(
  (dfa7$ID == 1 & ((dfa7$lower < 0.5 & dfa7$upper < 0.5)|(dfa7$lower > 0.5 & dfa7$upper > 0.5)))|
    (dfa7$ID != 1 & ((dfa7$lower < 0 & dfa7$upper < 0)|(dfa7$lower > 0 & dfa7$upper > 0))),"not","yes")

dfa7$ID <- factor(dfa7$ID,levels=c(1,2,3),labels=c("Marginal Means","AMCE","Diff. Marginal Means"))

dfa7$level <- fct_relevel(dfa7$level,"Much more strict",after=20)
dfa7$level <- fct_relevel(dfa7$level,"Constant",after=20)

dfa7$tech7 <- as.factor(as.character(dfa7$tech7))
dfa7$tech7 <- fct_relevel(dfa7$tech7,"High",after=2)

##Plot Results by Technology solves Issues

figa7 <- ggplot(dfa7, aes(x = level, y = predicted))+
  geom_hline(aes(yintercept=ifelse(ID=="Marginal Means",0.5,0)),size=1,color="grey")+
  geom_pointrange(aes(ymin=lower,
                      ymax=upper,
                      color=tech7,
                      shape=tech7,
                      alpha=sig),
                  position=position_dodge2(width=.75),
                  size=1.25,fatten=1)+
  facet_grid(feature~ID,
             scales="free",
             switch="y",
             space = "free")+
  scale_x_discrete(position = "top",
                   expand=c(-.65,0))+
  coord_flip()+
  theme_classic(base_family = "sans")+
  theme(text=element_text(size=16,family="Roboto Condensed"),
        strip.text.y.left=element_text(angle=0,hjust=0,
                                       vjust=0),
        strip.text.x=element_text(face="bold"),
        strip.placement.y="outside",
        axis.text.y = element_text(hjust=-1),
        panel.spacing.x=unit(1,"lines"),
        panel.background = element_rect(fill = NA, color = "grey"),
        legend.position="bottom",
        panel.grid.major = element_line(),
        legend.title=element_text(size=12))+
  labs(x="",y="")+
  scale_color_manual(values=c("#001524","#219ebc","#FF7D00"))+
  scale_alpha_manual(values=c(1,.25))+
  scale_shape_manual(values=c(rep(15,9)))+
  guides(alpha="none",color=guide_legend(title="tech7"),
         shape=guide_legend(title="tech7"))
figa7

ggsave(figa7,filename="Appendix_Figure7.PDF",width=12,height=7,dpi=500)


###############################################
##### Appendix Figure 8: Algorithm Anxiety ####
###############################################

CJa <- CJ[!is.na(CJ$ADMuneasy),]
CJa$ADMuneasy <- 6-CJa$ADMuneasy
CJa$ADMuneasy <- factor(CJa$ADMuneasy,levels=c(1,2,3,4,5),labels=c("Low","Low","Low","High","High"))


## Estimate Regression Model with interaction by ADMuneasy
ma8 <- setNames(lapply(att,function(x){
  svyglm(as.formula(paste("selected~",x,"*ADMuneasy")),design=svydesign(data=CJa,id=~ResponseId))
}),att)
## Estimate Marginal Means
mma8 <-  data.frame(bind_rows(setNames(lapply(att,function(x){
  ggeffect(ma8[[x]],terms=c(paste(x),"ADMuneasy"))
}),att),.id="feature"))%>%
  mutate(group = str_replace(group,"ADMuneasy",""))%>%
  rename(level=x,lower=conf.low,upper=conf.high,ADMuneasy=group)
## Estimate Average Marginal (Component) Effect by ADMuneasy
amcea8 <- bind_rows(setNames(lapply(att,function(x){
  summary(margins::margins(ma8[[x]],variables=paste(x),at=list(ADMuneasy=levels(CJa$ADMuneasy))))
}),att),.id="feature") %>% rename(level=factor,predicted=AME,ADMuneasy=ADMuneasy) %>%
  mutate(.,level=str_replace(level,feature,""))
##Create Dummy values for reference categories
b <- bind_rows(lapply(att,function(x){
  expand.grid(colnames(CJa[paste(x)]), levels(CJa[[paste(x)]]),levels(CJa$ADMuneasy))
}))
names(b) <- c("feature","level","ADMuneasy")
amcea8 <- merge(b, amcea8, by=c("feature","level","ADMuneasy"), all = T)
amcea8[is.na(amcea8)] <- 0
## Estimate Difference in Marginal Means between ADMuneasy
mmdiffa8 <-  data.frame(bind_rows(setNames(lapply(att,function(x){
  summary(margins::margins(ma8[[x]],variables=c("ADMuneasy"),at=setNames(list(lvls[[x]]), paste(x))))
}),att),.id="feature")) %>%
  ungroup() %>%
  mutate_at(att,as.character) %>%
  mutate(level = as.factor(as.character(coalesce(!!! rlang::syms(att)))),
         factor = str_replace(factor,"ADMuneasy","")) %>%
  select(.,-c(all_of(att))) %>% rename(predicted=AME,ADMuneasy=factor)
##Prepare for plotting
dfa8 <- bind_rows(mma8,amcea8,mmdiffa8,.id="ID")
levels(dfa8$level) <- lvl
dfa8$feature <- factor(dfa8$feature,levels=att,labels=ats)

dfa8$sig <- ifelse(
  (dfa8$ID == 1 & ((dfa8$lower < 0.5 & dfa8$upper < 0.5)|(dfa8$lower > 0.5 & dfa8$upper > 0.5)))|
    (dfa8$ID != 1 & ((dfa8$lower < 0 & dfa8$upper < 0)|(dfa8$lower > 0 & dfa8$upper > 0))),"not","yes")

dfa8$ID <- factor(dfa8$ID,levels=c(1,2,3),labels=c("Marginal Means","AMCE","Diff. Marginal Means"))

dfa8$level <- fct_relevel(dfa8$level,"Much more strict",after=20)
dfa8$level <- fct_relevel(dfa8$level,"Constant",after=20)

dfa8$ADMuneasy <- as.factor(as.character(dfa8$ADMuneasy))
dfa8$ADMuneasy <- fct_relevel(dfa8$ADMuneasy,"High",after=2)

##Plot Results by Algorithm Anxiety

figa8 <- ggplot(dfa8, aes(x = level, y = predicted))+
  geom_hline(aes(yintercept=ifelse(ID=="Marginal Means",0.5,0)),size=1,color="grey")+
  geom_pointrange(aes(ymin=lower,
                      ymax=upper,
                      color=ADMuneasy,
                      shape=ADMuneasy,
                      alpha=sig),
                  position=position_dodge2(width=.75),
                  size=1.25,fatten=1)+
  facet_grid(feature~ID,
             scales="free",
             switch="y",
             space = "free")+
  scale_x_discrete(position = "top",
                   expand=c(-.65,0))+
  coord_flip()+
  theme_classic(base_family = "sans")+
  theme(text=element_text(size=16,family="Roboto Condensed"),
        strip.text.y.left=element_text(angle=0,hjust=0,
                                       vjust=0),
        strip.text.x=element_text(face="bold"),
        strip.placement.y="outside",
        axis.text.y = element_text(hjust=-1),
        panel.spacing.x=unit(1,"lines"),
        panel.background = element_rect(fill = NA, color = "grey"),
        legend.position="bottom",
        panel.grid.major = element_line(),
        legend.title=element_text(size=12))+
  labs(x="",y="")+
  scale_color_manual(values=c("#001524","#219ebc","#FF7D00"))+
  scale_alpha_manual(values=c(1,.25))+
  scale_shape_manual(values=c(rep(15,9)))+
  guides(alpha="none",color=guide_legend(title="ADMuneasy"),
         shape=guide_legend(title="ADMuneasy"))
figa8

ggsave(figa8,filename="Appendix_Figure8.PDF",width=12,height=7,dpi=500)



############################################################
#### Appendix Figure 9: Estimate Attribute Interaction #####
############################################################

## Liability X Ethical Prioritization
ma9a <- svyglm(selected ~ Safety.Standards*Safety.Approval,design=svydesign(data=CJ,id=~ResponseId))
## Estimate Average Marginal (Component) Interaction Effect by Safety.Approval
amcea9a <- summary(margins::margins(ma9a,variables="Safety.Standards",
                                    at=list(Safety.Approval=levels(CJ$Safety.Approval))))%>%
  rename(level=factor,predicted=AME,group=Safety.Approval) %>%
  mutate(.,feature="Safety.Standards",level=str_replace(level,feature,""))
##Create Dummy values for reference categories
b <- expand.grid(colnames(CJ["Safety.Standards"]), levels(CJ[["Safety.Standards"]]),levels(CJ$Safety.Approval))
names(b) <- c("feature","level","group")
amcea9a <- merge(b, amcea9a, by=c("feature","level","group"), all = T)
amcea9a[is.na(amcea9a)] <- 0
amcea9a <- amcea9a %>% mutate(
  ID = "Safety Standards X\nSafety Approval By",
  feature="Safety Standards\nvs.current")


## Liability X Ethical Prioritization
ma9b <- svyglm(selected ~ Safety.Standards*Supervision,design=svydesign(data=CJ,id=~ResponseId))
## Estimate Average Marginal (Component) Interaction Effect by Safety.Approval
amcea9b <- summary(margins::margins(ma9a,variables="Safety.Standards",
                                    at=list(Supervision=levels(CJ$Supervision))))%>%
  rename(level=factor,predicted=AME,group=Supervision) %>%
  mutate(.,feature="Safety.Standards",level=str_replace(level,feature,""))
##Create Dummy values for reference categories
b <- expand.grid(colnames(CJ["Safety.Standards"]), levels(CJ[["Safety.Standards"]]),levels(CJ$Supervision))
names(b) <- c("feature","level","group")
amcea9b <- merge(b, amcea9b, by=c("feature","level","group"), all = T)
amcea9b[is.na(amcea9b)] <- 0
amcea9b <- amcea9b %>% mutate(
  ID = "Safety Standards X\nSupervision",
  feature="Safety Standards\nvs.current")


## Liability X Ethical Prioritization
ma9c <- svyglm(selected ~ Legal*Moral,design=svydesign(data=CJ,id=~ResponseId))
## Estimate Average Marginal (Component) Interaction Effect by Safety.Approval
amcea9c <- summary(margins::margins(ma9c,variables="Legal",
                                    at=list(Moral=levels(CJ$Moral))))%>%
  rename(level=factor,predicted=AME,group=Moral) %>%
  mutate(.,feature="Legal",level=str_replace(level,feature,""))
##Create Dummy values for reference categories
b <- expand.grid(colnames(CJ["Legal"]), levels(CJ[["Legal"]]),levels(CJ$Moral))
names(b) <- c("feature","level","group")
amcea9c <- merge(b, amcea9c, by=c("feature","level","group"), all = T)
amcea9c[is.na(amcea9c)] <- 0
amcea9c <- amcea9c %>% mutate(
  ID = "Liability X\nEthical Prioritization",feature="Liability for\nAccidents")

## Liability X Supervision
ma9d <- svyglm(selected ~ Legal*Supervision,design=svydesign(data=CJ,id=~ResponseId))
## Estimate Average Marginal (Component) Interaction Effect by Safety.Approval
amcea9d <- summary(margins::margins(ma9d,variables="Legal",
                                    at=list(Supervision=levels(CJ$Supervision))))%>%
  rename(level=factor,predicted=AME,group=Supervision) %>%
  mutate(.,feature="Legal",level=str_replace(level,feature,""))
##Create Dummy values for reference categories
b <- expand.grid(colnames(CJ["Legal"]), levels(CJ[["Legal"]]),levels(CJ$Supervision))
names(b) <- c("feature","level","group")
amcea9d <- merge(b, amcea9d, by=c("feature","level","group"), all = T)
amcea9d[is.na(amcea9d)] <- 0
amcea9d <- amcea9d %>% mutate(
  ID = "Liability X\nSupervision",
  feature="Liability for\nAccidents")


## Privacy X Ethical Prioritization
ma9e <- svyglm(selected ~ Data*Moral,design=svydesign(data=CJ,id=~ResponseId))
## Estimate Average Marginal (Component) Interaction Effect by Safety.Approval
amcea9e <- summary(margins::margins(ma9e,variables="Data",
                                    at=list(Moral=levels(CJ$Moral))))%>%
  rename(level=factor,predicted=AME,group=Moral) %>%
  mutate(.,feature="Data",level=str_replace(level,feature,""))
##Create Dummy values for reference categories
b <- expand.grid(colnames(CJ["Data"]), levels(CJ[["Data"]]),levels(CJ$Moral))
names(b) <- c("feature","level","group")
amcea9e <- merge(b, amcea9e, by=c("feature","level","group"), all = T)
amcea9e[is.na(amcea9e)] <- 0
amcea9e <- amcea9e %>% mutate(
  ID = "Access to Telemtry Data X\nEthical Prioritization",feature="Access to\nTelemtry Data")

## Privacy X Supervision
ma9f <- svyglm(selected ~ Data*Supervision,design=svydesign(data=CJ,id=~ResponseId))
## Estimate Average Marginal (Component) Interaction Effect by Safety.Approval
amcea9f <- summary(margins::margins(ma9f,variables="Data",
                                    at=list(Supervision=levels(CJ$Supervision))))%>%
  rename(level=factor,predicted=AME,group=Supervision) %>%
  mutate(.,feature="Data",level=str_replace(level,feature,""))
##Create Dummy values for reference categories
b <- expand.grid(colnames(CJ["Data"]), levels(CJ[["Data"]]),levels(CJ$Supervision))
names(b) <- c("feature","level","group")
amcea9f <- merge(b, amcea9f, by=c("feature","level","group"), all = T)
amcea9f[is.na(amcea9f)] <- 0
amcea9f <- amcea9f %>% mutate(
  ID = "Access to Telemtry Data X\nSupervision",
  feature="Access to\nTelemtry Data")

## Ethical Prioritization X Supervision
ma9g <- svyglm(selected ~ Moral*Supervision,design=svydesign(data=CJ,id=~ResponseId))
## Estimate Average Marginal (Component) Interaction Effect by Safety.Approval
amcea9g <- summary(margins::margins(ma9g,variables="Moral",
                                    at=list(Supervision=levels(CJ$Supervision))))%>%
  rename(level=factor,predicted=AME,group=Supervision) %>%
  mutate(.,feature="Moral",level=str_replace(level,feature,""))
##Create Dummy values for reference categories
b <- expand.grid(colnames(CJ["Moral"]), levels(CJ[["Moral"]]),levels(CJ$Supervision))
names(b) <- c("feature","level","group")
amcea9g <- merge(b, amcea9g, by=c("feature","level","group"), all = T)
amcea9g[is.na(amcea9g)] <- 0
amcea9g <- amcea9g %>% mutate(
  ID = "Ethical Prioritization X\nSupervision",
  feature="Ethical\nPrioritization")

plots <- lapply(list(amcea9a,amcea9b,amcea9c,amcea9d,amcea9e,amcea9f,amcea9g),function(df){
  ggplot(df,aes(x=level,y=predicted))+
    geom_hline(yintercept=0)+
    geom_pointrange(aes(ymin=lower,
                        ymax=upper,
                        color=group,
                        shape=group),
                    position=position_dodge2(width=.75),
                    size=1.25,fatten=7)+
    facet_grid(feature~ID,
               scales="free",
               switch="y",
               space = "free")+
    scale_x_discrete(position = "top",
                     expand=c(-.65,0))+
    coord_flip()+
    theme_classic(base_family = "sans")+
    theme(text=element_text(size=16,family="Roboto Condensed"),
          strip.text.y.left=element_text(angle=0,hjust=0,
                                         vjust=0),
          strip.text.x=element_text(face="bold"),
          strip.placement.y="outside",
          axis.text.y = element_text(hjust=-1),
          panel.spacing.x=unit(1,"lines"),
          panel.background = element_rect(fill = NA, color = "grey"),
          legend.position="bottom",
          panel.grid.major = element_line(),
          legend.title=element_blank())+
    labs(x="",y="")+
    scale_color_manual(values=c("#001524","#21a1bc","#FF7D00"))+
    scale_alpha_manual(values=c(1,.25))+
    scale_shape_manual(values=c("\u25BC","\u25C6","\u25B2","\u25A0"))+
    scale_y_continuous(breaks=seq(0,1,by=0.03))
})


figa9 <- cowplot::plot_grid(plots[[1]],plots[[2]],plots[[3]],plots[[4]],plots[[5]],plots[[6]],plots[[7]],ncol=2)
figa9
ggsave(figa9,filename = "Appendix_Figure9.PDF",height=14,width=10,dpi=500)


#####################################################################
#### Appendix Figure 10: Estimate Results by Climate vs Economy #####
#####################################################################

CJa <- CJ[!is.na(CJ$climateeco_bin),]

## Estimate Regression Model with interaction by climateeco_bin
ma10 <- setNames(lapply(att,function(x){
  svyglm(as.formula(paste("selected~",x,"*climateeco_bin")),design=svydesign(data=CJa,id=~ResponseId))
}),att)

## Estimate Marginal Means
mma10 <-  data.frame(bind_rows(setNames(lapply(att,function(x){
  ggeffect(ma10[[x]],terms=c(paste(x),"climateeco_bin"))
}),att),.id="feature"))%>%
  mutate(group = str_replace(group,"climateeco_bin",""))%>%
  rename(level=x,lower=conf.low,upper=conf.high,climateeco_bin=group)

## Estimate Average Marginal (Component) Effect by climateeco_bin
amcea10 <- bind_rows(setNames(lapply(att,function(x){
  summary(margins::margins(ma10[[x]],variables=paste(x),at=list(climateeco_bin=levels(CJa$climateeco_bin))))
}),att),.id="feature") %>% rename(level=factor,predicted=AME,climateeco_bin=climateeco_bin) %>%
  mutate(.,level=str_replace(level,feature,""))

##Create Dummy values for reference categories
b <- bind_rows(lapply(att,function(x){
  expand.grid(colnames(CJa[paste(x)]), levels(CJa[[paste(x)]]),levels(CJa$climateeco_bin))
}))
names(b) <- c("feature","level","climateeco_bin")
amcea10 <- merge(b, amcea10, by=c("feature","level","climateeco_bin"), all = T)
amcea10[is.na(amcea10)] <- 0

## Estimate Difference in Marginal Means between climateeco_bin
mmdiffa10 <-  data.frame(bind_rows(setNames(lapply(att,function(x){
  summary(margins::margins(ma10[[x]],variables=c("climateeco_bin"),at=setNames(list(lvls[[x]]), paste(x))))
}),att),.id="feature")) %>%
  ungroup() %>%
  mutate_at(att,as.character) %>%
  mutate(level = as.factor(as.character(coalesce(!!! rlang::syms(att)))),
         factor = str_replace(factor,"climateeco_bin","")) %>%
  select(.,-c(all_of(att))) %>% rename(predicted=AME,climateeco_bin=factor)


##Prepare for plotting
dfa10 <- bind_rows(mma10,amcea10,mmdiffa10,.id="ID")
levels(dfa10$level) <- lvl
dfa10$feature <- factor(dfa10$feature,levels=att,labels=ats)

dfa10$sig <- ifelse(
  (dfa10$ID == 1 & ((dfa10$lower < 0.5 & dfa10$upper < 0.5)|(dfa10$lower > 0.5 & dfa10$upper > 0.5)))|
    (dfa10$ID != 1 & ((dfa10$lower < 0 & dfa10$upper < 0)|(dfa10$lower > 0 & dfa10$upper > 0))),"not","yes")

dfa10$ID <- factor(dfa10$ID,levels=c(1,2,3),labels=c("Marginal Means","AMCE","Diff. Marginal Means"))

dfa10$level <- fct_relevel(dfa10$level,"Much more strict",after=20)
dfa10$level <- fct_relevel(dfa10$level,"Constant",after=20)

dfa10$climateeco_bin <- as.factor(as.character(dfa10$climateeco_bin))
dfa10$climateeco_bin <- fct_relevel(dfa10$climateeco_bin,"High",after=2)
##Plot Results by Sample

figa10 <- ggplot(dfa10, aes(x = level, y = predicted))+
  geom_hline(aes(yintercept=ifelse(ID=="Marginal Means",0.5,0)),size=1,color="grey")+
  geom_pointrange(aes(ymin=lower,
                      ymax=upper,
                      color=climateeco_bin,
                      shape=climateeco_bin,
                      alpha=sig),
                  position=position_dodge2(width=.75),
                  size=1.25,fatten=7)+
  facet_grid(feature~ID,
             scales="free",
             switch="y",
             space = "free")+
  scale_x_discrete(position = "top",
                   expand=c(-.65,0))+
  coord_flip()+
  theme_classic(base_family = "sans")+
  theme(text=element_text(size=16,family="Roboto Condensed"),
        strip.text.y.left=element_text(angle=0,hjust=0,
                                       vjust=0),
        strip.text.x=element_text(face="bold"),
        strip.placement.y="outside",
        axis.text.y = element_text(hjust=-1),
        panel.spacing.x=unit(1,"lines"),
        panel.background = element_rect(fill = NA, color = "grey"),
        legend.position="bottom",
        panel.grid.major = element_line(),
        legend.title=element_text(size=12))+
  labs(x="",y="")+
  scale_color_manual(values=c("#001524","#219ebc","#FF7D00"))+
  scale_alpha_manual(values=c(1,.25))+
  scale_shape_manual(values=c("\u25BC","\u25C6","\u25B2","\u25A0"))+
  guides(alpha="none",color=guide_legend(title="climateeco_bin"),
         shape=guide_legend(title="climateeco_bin"))
figa10

ggsave(figa10,filename="Appendix_Figure10.PDF",width=12,height=9,dpi=500)


#################################################################
### Appendix Figure 11:  Results by Modal Choice for Commute ####
#################################################################

CJa <- CJ[!is.na(CJ$commode),]

#Reduce Groups

CJa$commode <- as.factor(CJa$commode)
CJa$commode2 <- as.factor(
  case_when(CJa$commode == "Car(driver)" ~ "Commute Driving",
            CJa$commode == "Walking" | CJa$commode == "Bicycle" ~ "Commute Walk/Bike"))
CJa <- CJa[!is.na(CJa$commode2),]

## Estimate Regression Model with interaction by Modal Choice
ma11 <- setNames(lapply(att,function(x){
  svyglm(as.formula(paste("selected~",x,"*commode2")),design=svydesign(data=CJa,id=~ResponseId))
}),att)
## Estimate Marginal Means
mma11 <-  data.frame(bind_rows(setNames(lapply(att,function(x){
  as.data.frame(ggeffect(ma11[[x]],terms=c(paste(x),"commode2")))
}),att),.id="feature"))%>%
  mutate(group = str_replace(group,"commode2",""))%>%
  rename(level=x,lower=conf.low,upper=conf.high,commode2=group)
## Estimate Average Marginal (Component) Effect by Modal Choice
amcea11 <- bind_rows(setNames(lapply(att,function(x){
  summary(margins::margins(ma11[[x]],variables=paste(x),at=list(commode2=levels(CJa$commode2))))
}),att),.id="feature") %>% rename(level=factor,predicted=AME,commode2=commode2) %>%
  mutate(.,level=str_replace(level,feature,""))
##Create Dummy values for reference categories
b <- bind_rows(lapply(att,function(x){
  expand.grid(colnames(CJa[paste(x)]), levels(CJa[[paste(x)]]),levels(CJa$commode2))
}))
names(b) <- c("feature","level","commode2")
amcea11 <- merge(b, amcea11, by=c("feature","level","commode2"), all = T)
amcea11[is.na(amcea11)] <- 0
## Estimate Difference in Marginal Means between Modal Choice
mmdiffa11 <-  data.frame(bind_rows(setNames(lapply(att,function(x){
  summary(margins::margins(ma11[[x]],variables=c("commode2"),at=setNames(list(lvls[[x]]), paste(x))))
}),att),.id="feature")) %>%
  ungroup() %>%
  mutate_at(att,as.character) %>%
  mutate(level = as.factor(as.character(coalesce(!!! rlang::syms(att)))),
         factor = str_replace(factor,"commode2","")) %>%
  select(.,-c(all_of(att))) %>% rename(predicted=AME,commode2=factor)
##Prepare for plotting
dfa11 <- bind_rows(mma11,amcea11,mmdiffa11,.id="ID")
levels(dfa11$level) <- lvl
dfa11$feature <- factor(dfa11$feature,levels=att,labels=ats)

dfa11$sig <- ifelse(
  (dfa11$ID == 1 & ((dfa11$lower < 0.5 & dfa11$upper < 0.5)|(dfa11$lower > 0.5 & dfa11$upper > 0.5)))|
    (dfa11$ID != 1 & ((dfa11$lower < 0 & dfa11$upper < 0)|(dfa11$lower > 0 & dfa11$upper > 0))),"not","yes")

dfa11$ID <- factor(dfa11$ID,levels=c(1,2,3),labels=c("Marginal Means","AMCE","Diff. Marginal Means"))

dfa11$level <- fct_relevel(dfa11$level,"Much more strict",after=20)
dfa11$level <- fct_relevel(dfa11$level,"Constant",after=20)

##Plot Results by Modal Choice

figa11 <- ggplot(dfa11, aes(x = level, y = predicted))+
  geom_hline(aes(yintercept=ifelse(ID=="Marginal Means",0.5,0)),size=1,color="grey")+
  geom_pointrange(aes(ymin=lower,
                      ymax=upper,
                      color=commode2,
                      shape=commode2,
                      alpha=sig),
                  position=position_dodge2(width=.75),
                  size=1.25,fatten=2)+
  facet_grid(feature~ID,
             scales="free",
             switch="y",
             space = "free")+
  scale_x_discrete(position = "top",
                   expand=c(-.65,0))+
  coord_flip()+
  theme_classic(base_family = "sans")+
  theme(text=element_text(size=16,family="Roboto Condensed"),
        strip.text.y.left=element_text(angle=0,hjust=0,
                                       vjust=0),
        strip.text.x=element_text(face="bold"),
        strip.placement.y="outside",
        axis.text.y = element_text(hjust=-1),
        panel.spacing.x=unit(1,"lines"),
        panel.background = element_rect(fill = NA, color = "grey"),
        legend.position="bottom",
        panel.grid.major = element_line(),
        legend.title=element_text(size=12))+
  labs(x="",y="")+
  scale_color_manual(values=c("#001524","#219ebc","#FF7D00"))+
  scale_alpha_manual(values=c(1,.25))+
  scale_shape_manual(values=c(rep(15,8)))+
  guides(alpha="none",color=guide_legend(title="Predominant Mode of\nTransport for Commute"),
         shape=guide_legend(title="Predominant Mode of\nTransport for Commute"))
figa11

ggsave(figa11,filename="Appendix_Figure11.PDF",width=12,height=7,dpi=500)


##############################
### Regression Tables ########
##############################

stargazer::stargazer(m1, type = "text", out="Baseline.txt",
                     covariate.labels = c("Much more strict","Stricter",
                                          "Public Authority","Independent Body",
                                          "Choice of Driver","Protect Others",
                                          "Manufacturer",
                                          "Public Authorities","Companies",
                                          "Constant Supervision","After Warning Signal"),
          dep.var.labels   = "Probability Proposal was selected",
          single.row=T,
          column.labels =  c("Safety Standards","Oversight","Ethical Prio.",
                          "Liability","Access Data","Human Supervision"))



stargazer::stargazer(m5, type = "text", out = "Sample.txt",
                     dep.var.labels   = "Probability Proposal was selected",
                     single.row=T,
                     column.labels =  c("Safety Standards","Oversight","Ethical Prio.",
                                        "Liability","Access Data","Human Supervision"),
                     covariate.labels = c("Much more strict","Stricter",
                                          "Public Authority","Independent Body",
                                          "Choice of Driver","Protect Others",
                                          "Manufacturer",
                                          "Public Authorities","Companies",
                                          "Constant Supervision","After Warning Signal",
                                          "Japan","United States",
                                          "Much more strict : JP","Stricter : JP",
                                          "Much more strict : US","Stricter : US",
                                          "Public Authority : JP","Independent Body : JP",
                                          "Public Authority : US","Independent Body : US",
                                          "Choice of Driver : JP","Protect Others : JP",
                                          "Choice of Driver : US","Protect Others : US",
                                          "Manufacturer : JP",
                                          "Manufacturer : US",
                                          "Public Authorities : JP","Companies : JP",
                                          "Public Authorities : US","Companies : US",
                                          "Constant Supervision : JP","After Warning Signal : JP",
                                          "Constant Supervision : US","After Warning Signal : US"))


stargazer::stargazer(m4, type = "text", out = "TechnologyOptimism.txt",
                     dep.var.labels   = "Probability Proposal was selected",
                     single.row=T,
                     column.labels =  c("Safety Standards","Oversight","Ethical Prio.",
                                        "Liability","Access Data","Human Supervision"),
                     covariate.labels = c("Much more strict","Stricter",
                                          "Public Authority","Independent Body",
                                          "Choice of Driver","Protect Others",
                                          "Manufacturer",
                                          "Public Authorities","Companies",
                                          "Constant Supervision","After Warning Signal",
                                          "Tech.Opt.Med","Tech.Opt.High",
                                          "Much more strict : Tech.Opt.Med","Stricter : Tech.Opt.Med",
                                          "Much more strict : Tech.Opt.High","Stricter : Tech.Opt.High",
                                          "Public Authority : Tech.Opt.Med","Independent Body : Tech.Opt.Med",
                                          "Public Authority : Tech.Opt.High","Independent Body : Tech.Opt.High",
                                          "Choice of Driver : Tech.Opt.Med","Protect Others : Tech.Opt.Med",
                                          "Choice of Driver : Tech.Opt.High","Protect Others : Tech.Opt.High",
                                          "Manufacturer : Tech.Opt.Med",
                                          "Manufacturer : Tech.Opt.High",
                                          "Public Authorities : Tech.Opt.Med","Companies : Tech.Opt.Med",
                                          "Public Authorities : Tech.Opt.High","Companies : Tech.Opt.High",
                                          "Constant Supervision : Tech.Opt.Med","After Warning Signal : Tech.Opt.Med",
                                          "Constant Supervision : Tech.Opt.High","After Warning Signal : Tech.Opt.High"))


stargazer::stargazer(m2, type = "text", out = "EconomicLiberalism.txt",
                     dep.var.labels   = "Probability Proposal was selected",
                     single.row=T,
                     column.labels =  c("Safety Standards","Oversight","Ethical Prio.",
                                        "Liability","Access Data","Human Supervision"),
                     covariate.labels = c("Much more strict","Stricter",
                                          "Public Authority","Independent Body",
                                          "Choice of Driver","Protect Others",
                                          "Manufacturer",
                                          "Public Authorities","Companies",
                                          "Constant Supervision","After Warning Signal",
                                          "E.Liberalism.Med","E.Liberalism.High",
                                          "Much more strict : E.Liberalism.Med","Stricter : E.Liberalism.Med",
                                          "Much more strict : E.Liberalism.High","Stricter : E.Liberalism.High",
                                          "Public Authority : E.Liberalism.Med","Independent Body : E.Liberalism.Med",
                                          "Public Authority : E.Liberalism.High","Independent Body : E.Liberalism.High",
                                          "Choice of Driver : E.Liberalism.Med","Protect Others : E.Liberalism.Med",
                                          "Choice of Driver : E.Liberalism.High","Protect Others : E.Liberalism.High",
                                          "Manufacturer : E.Liberalism.Med",
                                          "Manufacturer : E.Liberalism.High",
                                          "Public Authorities : E.Liberalism.Med","Companies : E.Liberalism.Med",
                                          "Public Authorities : E.Liberalism.High","Companies : E.Liberalism.High",
                                          "Constant Supervision : E.Liberalism.Med","After Warning Signal : E.Liberalism.Med",
                                          "Constant Supervision : E.Liberalism.High","After Warning Signal : E.Liberalism.High"))


stargazer::stargazer(m3, type = "text", out = "CulturalConservatism.txt",
                     dep.var.labels   = "Probability Proposal was selected",
                     single.row=T,
                     column.labels =  c("Safety Standards","Oversight","Ethical Prio.",
                                        "Liability","Access Data","Human Supervision"),
                     covariate.labels = c("Much more strict","Stricter",
                                          "Public Authority","Independent Body",
                                          "Choice of Driver","Protect Others",
                                          "Manufacturer",
                                          "Public Authorities","Companies",
                                          "Constant Supervision","After Warning Signal",
                                          "Conservatism.Med","Conservatism.High",
                                          "Much more strict : Conservatism.Med","Stricter : Conservatism.Med",
                                          "Much more strict : Conservatism.High","Stricter : Conservatism.High",
                                          "Public Authority : Conservatism.Med","Independent Body : Conservatism.Med",
                                          "Public Authority : Conservatism.High","Independent Body : Conservatism.High",
                                          "Choice of Driver : Conservatism.Med","Protect Others : Conservatism.Med",
                                          "Choice of Driver : Conservatism.High","Protect Others : Conservatism.High",
                                          "Manufacturer : Conservatism.Med",
                                          "Manufacturer : Conservatism.High",
                                          "Public Authorities : Conservatism.Med","Companies : Conservatism.Med",
                                          "Public Authorities : Conservatism.High","Companies : Conservatism.High",
                                          "Constant Supervision : Conservatism.Med","After Warning Signal : Conservatism.Med",
                                          "Constant Supervision : Conservatism.High","After Warning Signal : Conservatism.High"))







#########################################################
### Results vs. Cregg for Differences between Samples ###
#########################################################

## Cregg

amcecregg <- cj(CJ,as.formula(paste("selected~",paste0(att,collapse="+"))),
                estimate="amce",
                id=~ResponseId,
                by= ~sample) %>% rename(.,predicted=estimate,sample = sample) %>% mutate(.,procedure="Cregg",
                                                                                         Analysis = "AMCE") %>%
  select(.,-c(1:3,7:9))

mmcregg <- cj(CJ,as.formula(paste("selected~",paste0(att,collapse="+"))),
              estimate="mm",
              id=~ResponseId,
              by= ~sample)%>% rename(.,predicted=estimate,sample = sample) %>% mutate(.,procedure="Cregg",
                                                                                      Analysis = "MM")%>%
  select(.,-c(1:3,7:9))

mmdiffcregg <- cj(CJ,as.formula(paste("selected~",paste0(att,collapse="+"))),
                             estimate="mm_diff",
                             id=~ResponseId,
                             by= ~sample)%>% rename(.,predicted=estimate,sample = sample) %>% mutate(.,procedure="Cregg",
                                                                                                     Analysis = "Diff.MM")%>%
  select(.,-c(1:3,7:9))

## Manual

amcemanual <- select(amce5,-c(5:7)) %>% mutate(.,procedure = "Manual",
                                               Analysis = "AMCE")
mmmanual <- select(mm5,-c(4)) %>% mutate(.,procedure = "Manual",
                                         Analysis = "MM")
mmdiffmanual <- select(mmdiff5,-c(4:6)) %>% mutate(.,procedure = "Manual",
                                                   Analysis = "Diff.MM")


df <- bind_rows(amcecregg,amcemanual,mmcregg,mmmanual,mmdiffcregg,mmdiffmanual,.id="ID")%>%select(.,-c("ID"))

levels(df$level) <- lvl
df$feature <- factor(df$feature,levels=att,labels=ats)

df_long <- df  %>%
  pivot_wider(names_from=c(procedure,Analysis),values_from=c(predicted,lower,upper)) %>%
  mutate(AMCE.Cregg.vs.Manual = round(predicted_Cregg_AMCE - predicted_Manual_AMCE,5),
         MM.Cregg.vs.Manual = round(predicted_Cregg_MM - predicted_Manual_MM,5)) %>%
  select(.,c(1:5,22,6:7,23))

df_long2 <- select(df_long,c(1:3,7:8))%>%
  group_by(sample,feature) %>%
  mutate(DiffMM_Cregg = predicted_Cregg_MM - predicted_Cregg_MM[1],
         DiffMM_Manual = predicted_Manual_MM - predicted_Manual_MM[1]) %>%
  ungroup()

df_long <- merge(df_long,df_long2) %>% mutate(Cregg.AMCE.vs.MM = predicted_Cregg_AMCE - DiffMM_Cregg,
                                              Manual.AMCE.vs.MM = predicted_Manual_AMCE - DiffMM_Manual)


write.csv(df_long, "est_comp.csv", row.names = F)
